#include "math.h"
#include "stdio.h"

int FFT_N = 8;

typedef struct float2{
    float r;
    float i;
    float2(float a=0, float b=0){
        r = a;
        i = b;
    }
}float2;

void swap(float2* arrayA, float2* arrayB, int size){
    float2 tmp;
    for (int i=0;i<size;i++){
        tmp = arrayA[i];
        arrayA[i] = arrayB[i];
        arrayB[i] = tmp;
    }
}

float2 operator +(float2 a, float2 b) 
{
    return float2(a.r+b.r, a.i+b.i);
}

float2 operator -(float2 a, float2 b) 
{
    return float2(a.r-b.r, a.i-b.i);
}

float2 operator *(float2 a, float2 b) 
{
    return float2(a.r*b.r-a.i*b.i, a.r*b.i+a.i*b.r);
}

void FFT2( float2* v ) {
    float2 v0 = v[0];
    v[0] = v0 + v[1];
    v[1] = v0 - v[1];
}

int expand(int idxL, int N1, int N2 ){
    return (idxL/N1)*N1*N2 + (idxL%N1);
}

void printArray(float2 *array, int size){
    for(int i=0;i<size;i++){
        printf("(%.1f,%.1f),", array[i].r, array[i].i);
    }
}

void FftIteration(int j, int N, int R, int Ns, float2* data0, float2*data1){
    float2 v[R];
    int idxS = j;
    float angle = -2*M_PI*(j%Ns)/(Ns*R);
    for( int r=0; r<R; r++ ) {
        v[r] = data0[idxS+r*N/R];
        float2 w = float2(cos(r*angle), sin(r*angle));
        v[r] = v[r]*w;
    }
    FFT2( v );
    int idxD = expand(j,Ns,R);
    for( int r=0; r<R; r++ ){
        data1[idxD+r*Ns] = v[r];
    }
}

float2* CPU_FFT(int N, int R, float2* data0, float2* data1) {
    for( int Ns=1; Ns<N; Ns*=R ) {
        for( int j=0; j<N/R; j++ ){
            FftIteration( j, N, R, Ns, data0, data1 );
            //printf("FftIteration Ns:%d j:%d\n", Ns, j);
            //printArray(data1, FFT_N);
        }
        swap( data0, data1 , FFT_N);
    }
    return data0;
}


int main(){
    printf("enter %s\n", __func__);
    float2 array0[FFT_N], array1[FFT_N];
    for(int i=0;i<FFT_N;i++){
        array0[i].r = i;
        array0[i].i = 0;
    }
    printArray(array0, FFT_N);
    CPU_FFT(FFT_N, 2, array0, array1);
    printf("after deal:\n");
    printArray(array0, FFT_N);
}